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Abstract. We present r-Java, an r-process code for open use, that performs r-process nucleosynthesis calculations. Equipped 
with a simple graphical user interface, r-Java is capable of carrying out nuclear statistical equilibrium (NSE) as well as static 
and dynamic r-process calculations for a wide range of input parameters. In this introductory paper, we present the motivation 
and details behind r-Java, and results from our static and dynamic simulations. Static simulations are explored for a range 
of neutron irradiation and temperatures. Dynamic simulations are studied with a parameterized expansion formula. Our code 
generates the resulting abundance pattern based on a general entropy expression that can be applied to degenerate as well as 
non-degenerate matter, allowing us to track the rapid density and temperature evolution of the ejecta during the initial stages 
of ejecta expansion. At present, our calculations are limited to the waiting-point approximation. We encourage the nuclear 
astrophysics community to provide feedback on the code and related documentation, which is available for download from the 
website of the Quark- Nova Project: |http://quarknova.ucalgary.ca7| 
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1. Introduction 

The rapid neutron capture process (r-process) is believed 
to be the mechanism for the nucleosynthesis of about half 



of the stable nuclei heavier than Iron (Burbidge et al. 1957 



ICameron 19571 ). Explosive and neutron-rich astrophysical en- 
vironments present ideal conditions for the r-process to take 
place. Possible candidate sites discussed in the literature 
include the neutrino-driven neutron-rich wind from proto- 
neutron stars ( Qian & Woosley 1996| 'Qian 2003), prompt ex- 
plosions of collapsed stellar cores, (Sumiyoshi et al. 2001, 



Wanajo et al. 2003 



neutron star decompression (Meyer 1989 



Goriely et al. 2004 1, tidal disruption in binary merger events 



IFreiburghaus et al. 1999a i, outflows in Gamma-Ray Bursts 
Surman et al. 2006l l etc. Most importantly, abundance data on 
r-process elements in metal-poor stars (Sneden et al. 2003) and 
certain radionuclides in meteorites (Qian & Wasserburg 2008) 
point toward the distinct possibility of multiple r-process 
sites, so that more than one of these sites may contribute 
to the observed abundance pattern of the r-process elements 
dTruran et al. 2002b . 



In 1985, an important work by Bethe & Wilson brought 
neutrino-driven winds from Type II Supernovae (SNe) to the 
forefront of the discussion about astrophysical sites for the r- 
process. Since then, much progress has been made in the mod- 
elling of type II SNe and neutrino winds of nascent neutron 
stars. Some examples of such work are Woosley et al. (1994), 
Takahashi et al. (1994), Qian & Woosley (1996), Cardall & 
Fuller (1997), Otsuki et al. (2000), Wanajo et al. (2001), 
Thomson et al. (2001). However, we have not so far arrived 
upon a single self-consistent astrophysical model that repro- 
duces the observed abundance of r-process elements and at the 
same time, is consistent with spectroscopic and meteoritic data. 
The most natural explanation is that the observed r-process pat- 
tern follows from a superposition of neutron capture events 
with differing neutron-to-seed ratios and exposure time-scales. 
Producing the third peak requires extreme values of entropy 
and dynamic time-scale that are a challenge for current models 
of Type II SNe explosions ( |Qian & Wasserburg 2008| l. Neutron 
star mergers provide a much larger neutron-to-seed ratio than 
type II SNe, but occur so infrequently that they fail to explain 
the early enrichment of r-process elements relative to Iron ob- 
served in metal-poor stars ( Qian 2000[ ). Clearly, further study 
is needed to better understand the r-process conditions, to this 
end we present a fairly general and flexible r-process code that 
is transparent and freely available to the nuclear physics and 
astrophysics community. 
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Fig. 1. Schematic representation of r- Java's capabilities and tools (see text for meaning of symbols). More details can be found 
in the manual available on the quark nova project website: |http : //quarknova.ucalgary. ca/| 



The motivation behind developing our own r-process code 
came in part from recognizing the difficulty entailed in having 
access to an openly shared cross-platform program where any 
interested researcher can examine and modify the dynamical 
evolution equations and the nuclear input. We perceive a dis- 
tinct need for this, as rare-isotope experiments such as FRIB 
later in this decade will make available new and improved nu- 
clear data for reaction rates relevant to the r-process. Besides, 
the r-process is of wide interest to the astrophysics community 
who might wish to explore different dynamical environments 
for the r-process, and not be as concerned, so to speak, about 
the input nuclear physics. To avoid a kind of "black-box" sit- 
uation arising from limited exchange of numerical codes, we 
decided to build a new code that is flexible enough to incorpo- 
rate different conditions and has detailed documentation. The 
code is capable of simulating nuclear statistical equilibrium 
(NSE), neutron bombardment of a static target as well as the 
expansion (including the relativistic case) of neutron-rich ma- 
terial. In this work, we discuss all of the above and also apply 
tfiis code to the unique environment of an ejected and expand- 
ing neutron star crust. For the code we designed a graphical 
user interface (GUI) which allows for real-time interpretation 
and analysis of simulations through data and graphs. We have 
made the r-process application freely available in Java-based 
format on the quark nova project websita4 where one can also 



' 'http://quarknova.ucalgary.ca/software/rJava/' . If using this code 
for calculations for a paper, we request that you please acknowledge 
this website and cite the online arXiv article, or the published version 
when it becomes available. 



find input nuclear data and supporting documentation on most 
aspects of our code. We request users to test it, and welcome 
questions and feedback. Our aim is to make our nucleosynthe- 
sis code, which we call r-Java, and its workings as transparent 
as possible so that any scientific user can adapt it to their needs. 

This paper is organized as follows. In section 2, we men- 
tion our sources for the nuclear input within r-Java and dis- 
play results for NSE that are in agreement with previous bench- 
marks. In sections 3 and 4, we present the simplified reaction 
network (waiting-point approximation) and simulation results 
for heavy-element production by static and dynamic neutron 
capture. In section 5, the conclusions from our simulations are 
discussed and we look ahead to future work. 

2. The R-process Code: r-Java 

Fundamental to our code is the waiting point approxima- 
tion which sets lower limits on the temperature of 2 x 10''K 
and neutron density of I0^° cm"^ for all our simulations 
dCameron et al. 1983 1). Each of the modules of r-Java; namely 
NSE, static and dynamic have their own appropriate set of in- 
put parameters and numerical solution methods (see Fig. [l]). 
For example, in the dynamic case, the user can choose the 
amount of heating from neutrinos, turn on/off nuclear fission, 
alpha or neutron decay and specify the density profile of the ex- 
panding material. The nuclear input data that is provided along 
with r-Java uses nuclear masses from experimental sources as 
much as possible (Audi & Wapstra 1995 ), and where none was 
available we used iMoller et al. 1997, For the two nuclear in- 
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puts: the neutron separation energy (S„) and the j0-decay rates 
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Fig. 2. Abundances of nucleons and selected nuclei in the NSE 
state with varying temperature, at p = 1 x 10^ g.cm^^ and = 
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Fig. 3. Abundances of selected nuclei in the NSE state with a 
varying electron-to-baryon fraction, at p = 1 x 10^ g.cm^^ and 
T = 3.5 
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Fig. 4. Abundances of nucleons and selected nuclei in the NSE 
state with varying density, at - 0.5 and T - 6.5 Tg 



{A^), we used ' Moller et al. 1997| A freedom provided to the 
users of r-Java is the ability to modify any nuclear input in or- 
der to simulate the effect that such a change would have on 
abundance yields. There are many different outputs that can 
be generated by r-Java as well as the ability to animate how a 
given parameter or abundance changes during a simulation run. 
A unique module included in r-Java displays the periodic table 
and in real-time shows which elements are being created and 
destroyed during a simulation. Also included is a module to 
study nucleosynthesis in a Quark-Nova dJaikumar et al. 2007] ), 
which will be discussed in detail in subsequent papers. The 
following subsections are dedicated to describing in detail the 
physics involved in each of the modules of r-Java and present- 
ing corresponding simulation results. 

2.1. Nuclear Statistical Equilibrium 

As a prelude to more complex situations, we begin by develop- 
ing a NSE code where the determining nuclear property is the 
binding energy. Consider an astrophysical plasma composed of 
photons, free neutrons and protons, and a mix of seed nuclei 
such as ^^Fe, with temperature and density high enough that 
the nuclear reactions assembling nuclei from free neutrons and 
protons are much faster than the expansion time scale. For large 
enough temperatures r9=r/(10'')K^ 6, the system would then 
reach NSE with forward and backward reaction rates equal. As 
long as such conditions prevail, the detailed balance equation 
holds 



(1) 



where yu, is the chemical potential of the species / (with 
Zi,Ni,Ai the corresponding atomic number, neutron number, 
and mass number, respectively). Assuming all the nuclei in the 
gas follow Maxwell-Boltzmann statistics, the particle number 
densities are 



(InkTnii 



(2) 



where B is the binding energy and g is the degeneracy fac- 
tor. In addition, mass and charge conservation imply that 



i 

Y^ZiYi^Y,. 



(3) 
(4) 



Combining Eqs. (j2]l and ([T]l with the conservation equations 
Q and Q, we obtain the following equations for the neutron 
and proton chemical potentials: 



2, A/n,(^p,/i„) -pA^A = 



(5) 



For any given Y^ (electron fraction), T (the temperature 
in Kelvin) and baryon mass density p (in g cm""*), we then 
solve this system by the Newton-Raphson method for //„ and 
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fip (Fig. [TJ, which we use to obtain the mass fractions of ev- 
ery nucleus through Eqn.([8]l, where Na is the Avogadro num- 
ber. The nuclei that dominate the abundances distribution in 
NSE depend only on p, T and B. High densities favor large 
nuclei due to the p-dependence while high temperatures fa- 
vor light nuclei since the Planck distribution contains ener- 
getic photons which photo-disintegrate heavy nuclei. We found 
the detailed balance method (Eqn. ([TJ) numerically amenable 
to the Newton-Raphson method. However, we should caution 
that the success of the Newton-Raphson method depends sen- 
sitively on a proper choice for the initial value of fip,fi„ and 
the starting composition. We chose to use the NSE calculations 
of 'Seite nzahl et al. 2008l as a test of our code. We have plotted 
our NSE results in Fig|2] We also display the mass fractions 
as a fu nction of Yg in Fig|3] Co mparing the former to the re- 
sults of Clifford & Tayler 1965 and the latter with Figure 3 of 
ISeiten zahl et al. 2008 we find very good agreement with these 
established results. In Fig|4] we have also shown a new plot: the 
mass fractions obtained with our code as a function of density, 
for fixed = 0.5 and J = 6.5 x lO'^K. 

3. Reaction Network and Static Case 

In general, a whole variety of reactions can occur simultane- 
ously, producing and destroying nucleus /, whose reaction in- 
duced density evolution is described by: 



(5) =i:«;o-Z«w. 

\ / p=const j j^i^ 



(6) 



where one-body reactions, such as decays and photo- 
disintegrations (eg. ry = ^i^j) and two body reactions {rjk - 
< (TV > jk rijiik where 6jk is the Kronecker delta) dom- 
inate. The individual N's in Eqn. (j6]l are positive or negative 
numbers specifying creation and annihilation. Then, the num- 
ber abundances F, and mass abundances X, are given by 



pNa 
niirii 



(7) 
(8) 



In terms of F,, the nuclear reaction network equation (|6]) 
becomes 



% - X N)AjYj + T^pNa < <Tv >j, YjY, 



(9) 



Applying this reaction network to the r-process, we will 
limit ourselves at the moment to the following interactions: 

- neutron capture 

- photo-disintegration 

- jS-decay 

For the static simulation, we assume a static, non- 
expanding chunk of neutron-rich material near drip densit 
(« 4 X 10"g/cc), described by the BPS equation of statf^ 



^ This leaves only T, n„ and Zq (the dominant nucleus at that den- 
sity) as inputs since is prescribed. 



(Baym et al. 1971b i, being bombarded by a large and constant 



flux of neutrons over a given period of time. Eqn. (j9]l is then 



Y(Z,A) = X'^_^^Y(Z-l,A) + <crv>n„Y{Z,A-\) 



+ Y(Z,A + 1)4^+1 - Y{Z,A){< o-v > tin 



(10) 



where {crv} is the thermally averaged neutron capture cross- 
section. Since we have a high neutron density (10^^-10^''cm"^) 
and high temperatures (r > 2 x lO^/T), neutron captures 
and photo-disintegrations occur much faster than j6-decays. For 
such conditions we can assume that (n,y) ^ (t, «) equilib- 
rium is established within every isotopic chain. We also assume 
steady flow while the high density and temperature conditions 



remain, viz.. 



F(Z,A) = 0, 



(11) 



This leads to steady y6-flow or the standard waiting-point ap- 
proximation in which we can express the abundance ratio of 
two neighboring isotopes within one isotopic chain in terms of 
only the neutron separation energy. In terms of the reduced re- 
action network ( fTO] ), we then have: 



y(z,A-H 1) _ 

Y(Z,A) " 

g(Z,A + l) IA + 1 
2g(Z,A) 



■ill 



2nh- 



\ljkTj 



,3/2 

e kT 



(12) 



Further, we introduce the abundance of an isotopic chain, Y{Z), 
defined as the sum of every isotope abundance in that chain: 



^'(•Z) = YjA Y(Z,A). With Eqn.( 12 1, we can now express every 
Y(Z,A) in terms of only the corresponding Y{Z) and popula- 
tion coefficients defined through Y(Z,A) = P{Z, A)Y(Z) which 
depend only on «„, T and S„. This allows us to dramatically 
reduce the number of equations from several thousand down 
to 136, the number of isotopic chains we have in our network. 
Now, we study the evolution of Y(Z) through y6-decay: 

Y(Z) = Yj^(Z,A) 

= Tj^4-i,a^(^- 1"^) - 4aY(Z,A)) . (13) 
This leads to 

Y(Z) = 2(4_; ^P(Z- I, A)) X Y(Z- 1) 

-2(4^^(2.^4)) xF(Z). (14) 

Introducing efl'ective j6-decay rates: A^^^ - ^^(/l^^f (Z, A)), we 
obtain a system of coupled differential equations: 



Y{Z) = 4{^F(Z - 1) - 4-^^F(Z) 



(15) 



As the abundances Y(Z) can vary by many orders of magni- 
tude within the integration interval, equation ( [T5] l is a very 
stiff system, so we use an implicit scheme, the Backward Euler 
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Fig. 5. Final abundances in the static simulations after one second of neutron bombardment and at different physical conditions. 
Top left: n„ = lO^^cm^^ and 7^9=2, Top right: n„ = lO^^^cm^^ and 7^9=2, Bottom left: n„ = lO^^'cm^^ and 7^9=4, Bottom right: 
«„ - 10^**cm"^ and Ti)-A. The solar abundance is shown in each case with vertical lines marking the location of the important 
peaks. 



Scheme, to treat the problem. To present the capabilities of r- 
Java in the static regime we varied the neutron density number 
and the temperature of the material to observe their impact on 
the r-process abundance (see Fig.jSj. 

Comparing the left and right figures in the top panel of 
Fig. |5] shows clearly that the third r-process peak at A=195 
is realized when the neutron number density is high (n„ - 
lO^^cm"^) and the temperature is close to 2 x lO^K. In real- 
ity, in dense neutron-rich environments neutron number densi- 
ties above lO^^cm"-* are possible, and we might expect large A 
nuclei to be copiously produced, but as we have not included 
fission in the static simulations, we have restrained our code at 
present to lower neutron densities (n„ - 10^^ - lO^^cm"-'), so 
that nuclei do not build up without limit. 

3.1. Static Simulation Results 

From the top panels of Fig. |5] it can be seen that our static 
simulations produce a weak third r-process peak at A=195 



even when the neutron number density is relatively high (n„ — 
lO^^cm"^). Higher neutron densities can produce a strong third 
peak (n„ ~ lO-'^cm"''), but the absence of fission presently lim- 
its our static code to moderate neutron densities. Absence of 
fission cycling also leads to a broad peak structure. The bot- 
tom panels of Fig. |5]show the effects of doubling the tempera- 
ture, relative to the top panel, keeping neutron density the same. 
As temperature increases heavy nuclei tend to disappear since 
photo-disintegrations of neutron-rich nuclei become highly ef- 
fective freezing the nuclear flow at lower A. In this case, the 
peaks in the mass distribution appear at lower mass numbers 
compared to the corresponding usual r-process peaks. This was 
found to be the case even for longer periods of neutron expo- 
sure, possibly suggesting that photo-disintegrations are domi- 
nating over neutron-captures within an isotopic chain, violat- 
ing the conditions of (n, y) equilibrium. Such a "shift" of the 
peaks to incorrect locations in mass number has been noticed in 
earlier dynamical situations, such as neutron star mergers (eg. 
|Freiburghaus et al. 1999b| l, where it is attributed to the choice 
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Fig. 6. Final abundances in the static simulations showing how the elements jS-decay back to stability once bombardment of 
neutrons stops. The cases shown are one second and one hundred seconds following the end of neutron bombardment. The left 
panel is the corresponding Z versus N distribution while the right panel displays how the abundances versus charge number 
evolve in time. 



of Yg. In our static simulations, does not cause the shift. 
Rather, in our case, this effect is driven by the high temperature, 
which favours photo-disintegrations and takes the r-process 
path closer to the valley of beta-stability, thereby downplay- 
ing the role of waiting-point nuclei in determining the abun- 



dance peaks (see Goriely & Arnould(1996) for the high neu- 
tron density and low temperature, Tg < 2, regime). Extremely 
high neutron densities are required to restore the peak to its 
correct position. 



3.2. The p-decay module 

Although we have implemented the waiting-point approxima- 
tion, r-Java also provides the option of evolving the elements by 
y6-decay back to stability from the time bombardment of neu- 
trons is abruptly terminated. Thus, one can isolate key y6-decay 
rates for particular reaction networks and studies of freeze-out, 
which is known to be important, for eg., in explaining the rare- 
earth peak near Europium ( ISurman et al. 1997] ). We use the 
Euler implicit scheme for each isotopic chain of elements with 
decay rate of each element exactly the ones we used to compute 
the effective rates defined in Eqn.( 15 i. Shown in Figure [6] are 
final abundances from the static simulation (under same input 
conditions as in Fig. [5]) y6-decaying one second and one hundred 
seconds after the bombardment of neutrons has stopped. An 
added convenience in r-Java is the fact that the y6-decay mod- 
ule can be run independent of nucleosynthesis calculations by 
specifying the initial composition as a vector. We added a mod- 
ification of the bi-diagonal solver, so that this solver accepts a 



vector of variable size. Currently the y8-decay module is limited 
to nuclei with mass between 50 and 300. 



4. Dynamic Case 

The dynamic module of r-Java allows the user to input the ini- 
tial temperature (To), density (po), electron fraction (Kg^o) and 
dominantly abundant element (Zo) of a chunk of material. The 
material then expands at a user specified expansion time-scale 
(t) for a given duration. In order for r-Java to accommodate a 
wide variety of ejection mechanisms we have allowed the den- 
sity profile (p (f)) to be freely modified by the use^ 

The relativistic expansion regime will be discussed in an 
upcoming paper, here we discuss the non-relativistic case. 
Users of r-Java are given the choice to include a simple fis- 
sion model for dynamic simulations, in which the last species 
of the network can be specified. For our dynamic simulations 
we assume that the last species in the network is Z-92, A-276 
which fissions into two fragments-one with Z=40 and one 
with Z-52. The code does not presently include other types 
of fission (neutron-induced, neutrino-induced, j6-delayed) or /3- 
delayed neutron emission. Further options given to the user of 
r-Java for dynamic simulations are the ability to specify the 



^ The default density profile in the dynamic module of r- 
Java comes from studies of non-relativistic decompression by 



[Meyer & Brown 1997| 



Pit) 



Po 



(1 + 



(16) 



Charignon et al.: r-process Nucleosynthesis with r-Java 



7 



fraction deposited from y6-decays, and allowing j6-decay and 
neutron decay to be turned on or off. 

Once the user provides the input parameters as above, the 
initial entropy/nucleon is determined (in ^^=1 units) from the 
initial density, electron fraction and temperature using 

Wn^T^ ulT 5 1 l gQ{2nmNkBT)^l^ \ 

s 1 1 1 In r , (17) 

45 p 3p 2)11^ \ n„h^ j 

where ^f-Qn^Y,, pY^^ is the chemical potential of the fully de- 
generate relativistic electrons, mM-939.\ MeV is the neutron 
mass and gQ-2 (spins) is the statistical weight of the neutron. 

Once the code begins to run, at each and every time step, 
we update the abundance of nucleus / if the fractional change 
dYil Yj < 0.1, else we decrease the time step until this condition 
is satisfied. From Eqn. ( 18 i, this implies that |5n,7n,- - dp/p\ < 
0.1. For non-relativistic expansion, both drii/rti and dp/p are 
small, but for relativistic expansion, dp/p can be large at early 
times. We can still ensure that the condition on the y,'s is sat- 
isfied by choosing our time step small enough that dp/p < 0.1, 
provided dp/p dominates over dni/ni. This is equivalent to an 
upper limit on the entropy which can be converted to an upper 
limit on the initial temperature rup(K)=0.01ny' Zo, where no 
is the initial ejecta number density (cm"-'). For example, tak- 
ing Zo = Zpi. (iron-rich ejecta), we find that an initial den- 
sity po ~ lO'-'g/cc requires that the initial temperature satisfy 
Tq < 10"K. If the temperature input by the user violates this 
density relation r-Java asks for a new choice of initial tempera- 
ture and shows the user the upper bound. Introducing heat into 
the system increases the entropy and can potentially reduce the 
maximum initial temperature for which our abundance crite- 
rion can be satisfied. However, heating due to nuclear transmu- 
tations can be ignored at early times in the case of fast relativis- 
tic expansion, since the expansion is much faster than typical 
beta-decay time-scales. Thus, the above estimate for the maxi- 
mum initial temperature remains valid for fast expansions. It is 
automatically satisfied for the non-relativistic case, where dp/p 
is much smaller. Furthermore, once satisfied at the outset, the 
abundance criterion dYi/Yi < 0.1 is seen to be always satisfied 
for all future times. These features are naturally built-in to our 
code, as detailed by the algorithm flow-chart seen in Figure |7] 

Equation ^ for the evolution of Y(Z) must now be mod- 
ified due to the density evolution, which was neglected before 
in the static calculations. The total derivative of F(Z) is 



Ye,o.Y{Zo)J 
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Fig. 7. Schematic representation of a single time step in our dy- 
namic code. In this chart po is the initial crust density, Yg Q is 
the initial electron to baryon ratio, Y(Zo) is the initial chemi- 
cal abundance of the crust. To is the initial temperature, Zo is 
the initial atomic number and no is the initial particle density 
of the crust. The neutron density (n„) is first found using the 
secant method and then reaction network is solved using the 
Crank-Nicholson method. After solving the reaction network 
evolution of the density (p{t)), heating due to /3 and neutron 
decay and nuclear fission is calculated. The neutron density is 
then re-evaluated and finally the maximal temperature change 
criterion is checked before moving on to the next time step. 



dYi 
dt 



= n; 



p 
dt 



1 / dui 



/ «,=const 



(18) 



p=const 



The first term, which describes changes in ejecta density, can be 
rewritten in term of 7, and then incorporated into our reaction 
network: 



d 



rij "p n,- 1 dp 
N^~dt ~ ~N^p^~dt 



1 dp 
"p~dt 



Yi. 



(19) 



For the dynamic simulation, we used a Crank-Nicholson 
method to solve Eqn.([T9]l. As in the static case, we were able to 



solve the linear system of abundance equations to find the vari- 
ations in abundances within a timestep. These nuclear transmu- 
tations generate entropy and lead to heating, so to keep track of 
this, we inverted the entropy from these changes to compute 
the change in temperature. This procedure was implemented in 
a self-consistent manner, i.e, the new temperature was fed back 
to check the value of entropy. The heating of the material via 
y6-decay is added to Eqn.( 17 1 as 6s - Sq/T and this whole ex- 
pression is inverted to get the new temperature after each time 
step. We check that the temperature remains above 2xlO^K at 
all times during r-process, so that the waiting point approxima- 
tion we employ is still valid. 
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Fig. 8. The final abundances (solid red line) for dynamic simulations (one second duration) of nucleosynthesis from expansion of 
a chunk of pure Iron situated in neutron-rich matter with Yg = 0.03. Solar abundances are shown as black crosses. For each panel 
the same initial temperature (4 x lO'^K) and density (1 x 10"g/cc) is used. The left column of panels use an expansion time-scale 
of T = 0.1s and for the right column of panels the expansion time-scale is t = 0.001s. Each row considers a different density 
profile: top p{t) oc {t/Ty\ middle p(t) oc (t/ry^, bottom p(f) oc {t/Ty^. 



To complete the time step, we also need to explicitly calcu- 
late the neutron density, which we can do using (n, y) ^ (y, n) 
equilibrium, as long as temperatures and densities remain high 
enough. This ensures equilibrium within any isotopic chain, 
which will instantly redistribute Y{Z) between all isotopes in 
accordance with the neutron separation energy. Once again, we 
use the population coefficients which depend on the neutron 
density, so to get the equilibrium value of neutron density at 



each timestep, we solve the implicit equation of mass conser- 
vation: 

n„m„ ^pNa{l-Yj ^(««)) ' (20) 

where 2 A is the mass of all nuclei except neutrons. This quan- 
tity depends on the neutron density through the population co- 
efficients, which give the effective mass of every isotopic chain. 
We used a secant method to solve this equation and get the 
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Fig. 9. The final Z versus N distribution is plotted for each dynamic simulation (after one second duration). In each panel the 
solid green line is at = Z, and the large crosses mark the location of the closed neutron shells. For each panel the same initial 
temperature (4 x lO^K), density (1 x 10''g/cc) and electron fraction (0.03) are used. The left column of panels use an expansion 
time-scale of t = 0.1s and for the right column of panels, the expansion time-scale is t = 0.001s. Each row considers a different 
density profile: top p(t) oc (t/Ty \ middle p(t) oc (f/r)"^, bottom p(t) oc (t/ry^. 



equilibrium value of n„. Then the system is completely de- 
scribed and the code can move to the next time step. We can run 
the code as long as the waiting point approximation is valid, for 
«„ > lO^^cm and T > 2x lO'^K (see algorithm in Figure[7]i. 



4.1. Dynamic Simulation Results 

In Fig. [8] we explore the effect on final abundance of changing 
the density profile of the expanding material. To do this we 
considered the following general density profile; 



(21) 
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Time (s) Time (s) 

Fig. 10. The evolution of temperature over the simulation run is plotted for each dynamic simulation (50% heating from j6-decays 
is assumed). For each panel the same initial temperature (4 x lO'^K), density (1 x 10"g/cc) and electron fraction (0.03) are used. 
The left column of panels use a expansion time-scale of t = 0.1s and for the right column of panels the expansion time-scale is 
T = 0.001s. Each row considers a different density profile: top p(f) oc (f/r)"', middle p(t) oc {t/Ty^, bottom p(t) oc (t/ry^. 



For Fig. [8] Fig.[9]and Fig. 10 in the top panels a - 1, the middle 
panels a = 2 and the bottom panels a - 3, where increasing 
values of a roughly simulate increasing degrees of freedom for 
expansion. The left column of these figures use an expansion 
time-scale of t = 0.001s while for the right column t - 0.1s. 
We left the remaining input parameters the same for all panels 
(To = 4 X lO^K, Ye_o=0.03, Zo=26 and a simulation duration of 
one second). From the left column of Fig.[8]it can be seen that 
the density drops off too rapidly to produce elements heavier 



than A=120 in all but the ff = 1 case. This is because the expan- 
sion is too fast for neutron captures to occur for a sufficiently 
long period of time to make heavy elements. An opposite trend 
can be seen in the t = 0.1s cases where the heavy element pro- 
duction increases with the degrees of freedom for expansion. 
This is because for slow expansions, the density, particularly 
for low degrees of freedom (a -I) remains high for an extended 
period of time, blocking beta-decays that are needed to move 
the flow up to higher mass numbers. Consequently, the third 



Charignon et al.: r-process Nucleosynthesis with r-Java 



11 



peak is absent or weak. For a-2 or 3, the density fall off is 
rapid enough to compensate this effect, so a strong 3rd peak 
is seen. The effect of the beta-decay blocking can be clearly 
seen in the Z vs. N plots shown in Fig.|9] In the top panel, right 
column of Fig. |9] it is clear that the slower expansion case al- 
lows elements to build up beyond the magic numbers, while for 
faster expansion (top-left panel of Fig.|9]) or multiple degrees of 
freedom (right column of 2nd and 3rd panel) there exists clear 
steps at each magic number, revealing the imprint of the extra 
stability. It is interesting to note that the final abundances seen 
in the top left (r = 0.001, a = 1) and bottom right (t = 0.1, 
a - J) panels of both Fig. [8] and |9] are remarkably similar. 
This similarity is due to the fact that after one second of ex- 
pansion the effects of the different a and r values will directly 
offset one another leaving the final density of the expanding 
chunk the same. The slight difference in final abundances of 
the two simulations are due to temperature effects. Finally, Fig. 
[To] displays the temperature evolution for each of our dynami- 
cal simulations. The left column of Fig. [TO] exemplifies how in 
the short expansion time-scale adiabatic cooling dominates the 
heating from beta-decays, while in the longer expansion time- 
scale cases (the right column of Fig. 10 1 beta decays have a 



strong heating effect early on, specially for a-2,3. 

5. Discussion and Conclusion 

We have constructed and tested a nucleosynthesis code for the 
r-process under static and dynamic conditions. Called r-Java, 
it is made available online through an easy-to-use interface at 
|http://quarknova.ucalgary.ca/] . The user is given a number of 
input options as well as the ability to easily modify the nuclear 
data used in the code. In this paper, we have studied NSE, as 
well as both static and dynamic simulations. 

Exploring static simulations for a range of neutron irra- 
diation and temperatures, we found that neutron density ~ 
lO^^cm"^ and above is required for a strong third r-process 
peak. With increasing temperature beyond 7^9=2, nuclei photo- 
disintegrate extremely fast, before significant pile-up can oc- 
cur at the true waiting-point nuclei. This resulted in abundance 
peaks at smaller mass numbers. This feature is compensated 
by increasing the neutron density; however, the effect of fis- 
sion will then play a role in determining the final abundances. 
We intend to include fission-rates in our static simulations in 
the future, so that higher neutron densities can be reliably ad- 
dressed. We have included a y6-decay module that allows the 
user to track the flow of y6-decays once neutron bombardment 
stops. This should be useful in studying freeze-out conditions 
and quasi-equilibrium reaction networks. 

Using the dynamic module of r-Java we compared differ- 
ent density profiles for two expansion time-scales. We found 
that at the shorter expansion time-scale (t = 0.001s) only the 
slowest evolving density profile (pit) oc poit/ry^) produced 
a significant abundance of elements above A=120. For the 
longer expansion time-scale (t - 0.1s) the heavy abundance 
yield increased with faster evolving density profiles, where the 
strongest A=195 peak occurred in the simulation that consid- 
ered pit) oc po (t/r)^^. At present, all our calculations are car- 
ried out in the waiting-point approximation. 



The large qualitative differences in the r-process pattern we 
observe as a function of dynamical expansion parameters are 
important for alternate sites of the r-process, such as neutron 
star mergers or spherically symmetric decompressing neutron 
star matter In the case of mergers, tidal forces can lead to a 
"tube of toothpaste" quasi-lD ejection, whereas spherical de- 
compressions resemble the 3D case. From Fig[8j we see that 
fast expansion is required if heavy elements are made in merg- 
ers, while less energy (slower expansion) is required for spher- 
ical decompression. This is the case, for eg., in Quark-Novae, 
which will be discussed in a follow-up paper We also note the 
curious fact that the peak at A = 3, corresponding to Tritium, 
is formed from the sequence of reactions comprised of neutron 
decay followed by formation of Deuterium and then Tritium. 
Since the half-life of Tritium is relatively short (12.32 years), a 
large amount of its daughter element, He-3, would be formed 
and dispelled into the surrounding space. In the case of a dual- 



shock Quark-Nova ( |Leahy & Ouyed 2008|[Ouyed et al. 20101 ), 
this can have important consequences, through spallation on 
supernova ejecta, for cosmic rays and anomalous abundance 
ratios in some supernovae (eg., Cas A). 

Our future efforts are aimed in three directions: (1) ap- 
plications of r-Java to Neutron Star Mergers and the Quark- 
Nova scenario. The Quark-Nova is an ejection mechanism for 
neutron star material powered by the phase conversion from 
hadronic to deconfined quark matter occurring at the core of 
a neutron star ( [Ouyed et al. 2002] IKeranen et al. 2003t . The 
extreme energetics and high neutron-to-seed ratio makes the 
Quark-Nova an ideal candidate for an astrophysical r-process 
site (IJaikumar et al. 2007l l: (2) including the effects of freeze- 
out through a set of quasi-equilibrium reaction networks and 



other processes neglected in Eq.( 13 1 such as beta-delayed neu- 
tron emission, to obtain a more realistic pattern of final abun- 
dances, which can then be used for schematic chemical evolu- 



tion studies; (3) based on Niebergal et al. 2010 , further numeri- 
cal study of the dynamics of the nuclear-quark conversion front 
inside the neutron star to better constrain the Lorentz factor of 
the ejecta, and possibly develop a semi-analytic approach to 
constrain the parameter space, similar to what has been done 
for neutrino-driven winds from supernovae. 
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